* ********************************************************************************************
* The Effects of Weather Shocks on Economic Activity: What are the Channels of Impact?
* Sebastian Acevedo, Mico Mrkaic, Natalija Novta, Evgenia Pugacheva, Petia Topalova
* With support from Gavin Asdorian, Olivia Ma, Jilun Xing and Yuan Zeng 
* Replication files for Table 4; Figure 8
* ********************************************************************************************
clear all
set more off
set matsize 11000

cd "C:\Users\..\replication" // set your working directory to where the replication folder is saved

* the Y variable that will be used in regressions
local Y_list ava_rlcu_wdi ag_prd_crop mva_rlcu_wdi sva_rlcu_wdi ni_r nm_r mrt_i hdi

* Horizons
global k 7

* ****************************************************************************************
* Data
* ****************************************************************************************
use "data/subnational_dataset.dta", clear

drop if year > 2016

* interaction with AE dummy for temperature and precipitation
foreach X of varlist pwtemp pwprecip {
	gen `X'_ae = ae * `X'
	label var `X'_ae "Interaction ae * `X'"
}

* temperature, keep sample within [15; 20]
bys provcode: egen mn_pwtemp = mean(pwtemp)
keep if inrange(mn_pwtemp, 15, 20)

foreach X in pwtemp pwprecip {
	gen `X'2=`X'^2
}

gen sqyear=year^2

by provcode (year), sort: gen d0_ln_gdppc = ln_gdppc - l.ln_gdppc
forval t = 1/$k {
	by provcode (year), sort: gen d`t'_ln_gdppc = f`t'.ln_gdppc - l.ln_gdppc
}

* Create WDI region year fixed effects
gen tempvar1 = wdi_region + " " + string(year) if wdi_region != ""
encode tempvar1, generate(ry)
drop tempvar1

xi i.ifscode*year i.ifscode*sqyear  i.ry i.year
* interact the RY with the AE dummy
xi i.ae*i.ry
gen l_d0_ln_gdppcl_ae = ae* l.d0_ln_gdppc

* ****************************************************************************************
* Variables where to store results
* ****************************************************************************************
gen period = -1 in 1
foreach X in pwtemp pwprecip pwtemp_ae pwprecip_ae lincom_ae{
	gen `X'_coef = 0 in 1
	gen `X'_se = 0 in 1
	gen `X'_pvalue = . in 1
}

gen nobs = . in 1
gen nprov_ae = . in 1
gen nprov_non = . in 1
gen nprov_total = . in 1
gen ncountry_ae = . in 1
gen ncountry_non = . in 1
gen ncountry_total = . in 1
gen year_start = . in 1
gen year_end = . in 1
gen r2 = . in 1

* ****************************************************************************************
* Program
* ****************************************************************************************
forval t = 0/$k {
	areg d`t'_ln_gdppc pwtemp pwprecip pwtemp_ae pwprecip_ae l.pwtemp l.pwprecip l.pwtemp_ae l.pwprecip_ae l_d0_ln_gdppcl_ae l.d0_ln_gdppc _Iae_* _Iry_* _IaeXry_*, cluster(provcode) absorb(provcode)
	
	foreach X in pwtemp pwprecip pwtemp_ae pwprecip_ae {
		replace `X'_coef = _b[`X'] in `=`t'+2'
		replace `X'_se = _se[`X'] in `=`t'+2'
		replace `X'_pvalue = 2 * ttail(e(df_r),abs(_b[`X']/_se[`X'])) in `=`t'+2'
	}

	replace period = `t' in `=`t'+2'

	replace nobs = `e(N)' in `=`t'+2'

	* number of provinces in AE
	levelsof provcode if e(sample) & ae == 1, local(province_list)
	local province_count: list sizeof province_list
	replace nprov_ae = `province_count' in `=`t'+2'

	* number of provinces in NON AE
	levelsof provcode if e(sample) & ae == 0, local(province_list)
	local province_count: list sizeof province_list
	replace nprov_non = `province_count' in `=`t'+2'

	replace nprov_total = `e(N_clust)' in `=`t'+2'

	* number of countries
	levelsof ifscode if e(sample), local(country_list)
	local country_count: list sizeof country_list
	replace ncountry_total = `country_count' in `=`t'+2'

	* number of countries in AE
	levelsof ifscode if e(sample)  & ae == 1, local(country_list)
	local country_count: list sizeof country_list
	replace ncountry_ae = `country_count' in `=`t'+2'

	* number of countries in NON AE
	levelsof ifscode if e(sample)  & ae == 0, local(country_list)
	local country_count: list sizeof country_list
	replace ncountry_non = `country_count' in `=`t'+2'

	* Years of the sample
	qui sum year if e(sample)
	replace year_start = r(min) in `=`t'+2'
	replace year_end = r(max) in `=`t'+2'

	* r squared
	replace r2 = e(r2) in `=`t'+2'

	lincom _b[pwtemp] + _b[pwtemp_ae]
	replace lincom_ae_coef = r(estimate) in `=`t'+2'
	replace lincom_ae_se = r(se) in `=`t'+2'
	replace lincom_ae_pvalue = 2 * ttail(r(df),abs(r(estimate)/r(se))) in `=`t'+2'
}

export excel period - r2 in 1/9 using "output/Table_4.xlsx", sheet("column_4") sheetreplace firstrow(var)
